/*******************************************************************************

	SSC0143 - Programacao Concorrente
	Prof. Julio Cesar Estrela
	
	T1 - Calculo de PI com 10000000 de casas decimais
	
	Grupo: 
	
	Daniel de Oliveira Pereira		6427171
	Guilherme Takasawa Yagui		6426698
	Vinicius Marques Stocco			6427233
	
	Março de 2013
	
********************************************************************************	
				Algoritmo de Borwein - Versao Sequencial	
*******************************************************************************/	
	
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <gmp.h>

// numero de casas decimais corretas de PI
#define PREC_DGT	10000000

// numero de bits minimo para armazenar PI com o numero de casas decimais corretas 
#define PREC_BIT 	(int) ceil((log(10) / log(2)) * PREC_DGT)

// numero de iteracoes necessario para calcular PI com o numero de casas decimais corretas desejado 
#define ITR 		(int) ceil(log(PREC_DGT) / log(4))

int main(char argc, char *argv[])
{
	FILE *f;
	
	clock_t t;					// tempo

	mpf_t ak, yk, ak1, yk1, bk; // calculo de pi
	mpf_t tmp1, tmp2;			// auxiliares
	
	unsigned int i;				
	unsigned int aux;
	
	// define a precisão desejada para cada variavel
	mpf_set_default_prec(PREC_BIT);
	
	// marca  o clock inicial da execucao do algoritmo
	t = clock();
	
	f = fopen("borwein_sequencial.txt", "w");
	
	// aloca memoria para as variaveis utilizadas para calcular PI
	mpf_init(ak);
	mpf_init(yk);
	mpf_init(ak1);
	mpf_init(yk1);
	mpf_init(bk);
	mpf_init(tmp1);
	mpf_init(tmp2);
	
	// inicializa a(0)
	mpf_sqrt_ui(tmp1, 2);
	mpf_mul_ui(tmp1, tmp1, 4);
	mpf_ui_sub(ak, 6, tmp1);
	
	// inicializa y(0)
	mpf_sqrt_ui(tmp2, 2);
	mpf_sub_ui(yk, tmp2, 1);
	
	
	for (i = 0; i < ITR; i++)
	{
		aux = (unsigned int) pow(2.0, (2.0 * i + 3.0));
		
		// calcula y(k+1)
		// y(k+1) = (1-(1-y(k)^4)^(1/4))/(1+(1-y(k)^4)^(1/4))
		mpf_pow_ui(tmp1, yk, 4);
		mpf_ui_sub(tmp1, 1, tmp1);
		mpf_sqrt(tmp1, tmp1);
		mpf_sqrt(tmp1, tmp1);
		mpf_set(tmp2, tmp1);
		mpf_ui_sub(tmp1, 1, tmp1);
		mpf_add_ui(tmp2, tmp2, 1);
		mpf_div(yk1, tmp1, tmp2);
		
		// calcula a(k+1)
		// a(k+1) = ak*(1+y(k+1))^4-2^(2k+3)*y(k+1)*(1+y(k+1)+y(k+1)^2)
		mpf_add_ui(tmp1, yk1, 1);
		mpf_pow_ui(tmp1, tmp1, 4);
		mpf_mul(tmp1, tmp1, ak);
		mpf_pow_ui(tmp2, yk1, 2);
		mpf_add(tmp2, tmp2, yk1);
		mpf_add_ui(tmp2, tmp2, 1);
		mpf_mul(tmp2, tmp2, yk1);
		mpf_mul_ui(tmp2, tmp2, aux);
		mpf_sub(ak1, tmp1, tmp2);
		
		// calcula PI = b(k) =  1/a(k+1)
		mpf_ui_div(bk, 1, ak1);
		
		//gmp_printf("\n %.*Ff", PREC_DGT, bk);
		gmp_fprintf(f, "%.*Ff \n", PREC_DGT, bk);
		
		mpf_set(yk, yk1);
		mpf_set(ak, ak1);
	}
	printf("\n Numero de iteracoes : %d \n", ITR);
	
	// libera a memoria alocada para calcular PI
	mpf_clear(ak);
	mpf_clear(yk);
	mpf_clear(ak1);
	mpf_clear(yk1);
	mpf_clear(bk);
	mpf_clear(tmp1);
	mpf_clear(tmp2);
	
	fclose(f);
	
	// calcula o numero clocks necessario para executar o algoritmo
	t = clock() - t; 
	
	printf("\n Tempo de execucao : %.3f s \n\n", (double) t / CLOCKS_PER_SEC);

	return 0;
}
